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ABSTRACT 

The incorporation of algebraic turbulence models in a solver for the 2D compressible 
Navier-Stokes equations using triangular grids is described. A practical way to use the 
Cebeci Smith model, and to modify it in separated regions, is proposed. The ability of the 
model to predict high speed, perfect gas boundary layers is investigated from a numerical 
point of view. 
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NOMENCLATURE 


Cr skin friction coefficient 

1 Pootiio 

C h heat flux coefficient -- H -- 

Pootioo 

C p pressure coefficent 2 

Z? deformation tensor £>(u) = Vu + Vu* — |V • u/ 
E total energy 

I identity operator 

k kinetic energy of turbulence 

M Mach number 

P pressure 

Pr Prandtl number ^ 

q heat flux 

Re free stream Reynolds number 

Reg momentum thickness Reynolds number e '^' 9 
T temperature 

u velocity vector u = (u, v ) 

u component of the velocity tangent to the wall 
u, speed scale 


in 


W vector of conserved variables 



x coordinate tangent to the wall 

y coordinate normal to the wall 

a Clauser’s constant a = 0.0168 

6 boundary layer thickness defined by u/u e = 0.99 

6* displacement thickness / 0 5 (1 — -f^)dy 

Si incompressible displacement thickness /g (1 — ^)dy 

e isotropic part of turbulence energy dissipation 

A heat conductivity 

H molecular viscosity 

k Von Karman’s constant k = 0.4 

"V kinematic viscosity ^ 

V gradient operator 

V • divergence operator 
<jj vorticity 
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SUBSCRIPTS 


b at backflow edge 

c at boundary layer edge 

i at streamwise position iAx 

j at crossflow position jAy 

t turbulent 
w at wall 

oo free stream conditions 

SUPERSCRIPTS 

' fluctuating quantities 
- averaged 
t transpose 
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I. INTRODUCTION 


Solving numerically the laminar compressible Navier-Stokes equations for practical con- 
figurations is now within reach [Ca], [Chi]. Flows around complete aircrafts, including an 
impressive high Mach number flow around a shuttle [Jl], have been successfully computed. 
However, the complexity of the configurations of industrial interest make these computa- 
tions to remain a challenging task, even when using modem supercomputers. The results 
obtained are for steady and laminar flows and still make full use of the computer capability 
available. This, together with the concern of affordability, is why cost is the major issue 
when more realistic physical models are taken into account. Still, it is well-known that 
most configurations of interest are at least locally turbulent, and that the modelling of 
this turbulence is critical for the reliability of the computations. Though a both cheap 
and reasonably accurate turbulent closure model is an unavoidable part of a compressible 
viscous code. 

The study of the second part of the reentry of the European shuttle Hermes prompts 
interest in high-speed, moderate Reynolds flows. With Mach numbers ranging from 5 to 
10, real gas effects can be neglected in a first approximation. The purpose of the study is 
not to predict the details of the flow, which are probably beyond reach anyway, but to give 
reasonable estimates (say ±10%) of the pressure, friction, and heat transfer, and through 
those of the flight performances of the shuttle. 

An implicit algorithm to solve the 2D (3D) compressible Euler equations on unstruc- 
tured triangular (tetrahedral) grids has been developed in the recent years [Stl, St2] . It 
is based on a finite volume approximation of the equations in conservation form, using 
control volumes defined by the medians of the triangle (tetrahedras) surrounding each 
node, and Osher’s flux formula [Osl]. It has been extended to solve the Reynolds aver- 
aged Navier-Stokes equations, and showed to be able to compute efficiently and accurately 
two-dimensional laminar boundary layers [Rol]. The purpose of this work is to include in 
it some relevant turbulence model. 

Our concern here will be two dimensional compressible boundary layers and compres- 
sion corners, which are the two relevant simple flows (see Fig. 1). Numerous turbulence 
models have been proposed to solve these problems, and a review would make a (tough) 
paper by itself. However, they can be classified by the number of extra P.D.E. they intro- 
duce, by the order of the closing relation, and by its linearity, i.e., by whether the model is 
of the eddy viscosity type or not. Nonlinear models have been shown to be an improvement 
over liner ones in a number of situations (see e.g., [Spl]), but they will not be our concern 
here, since we restrict ourselves to two dimensional situations in this first approach. 

The choice is then on the number of PDE one wishes to solve; it is our opinion that 
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1 and 2 equation models offer no significant improvement over algebraic models for the 
computations of the near-wall flows which are our concern. Consequently, we will restrict 
ourselves to algebraic models, in which we may later introduce some streamwise history 
effect through an explicit relaxation model, patterned after that of Shang and Hankey [Shi]. 
For attached flows, the models of Cebeci and Smith [Cel] and of Baldwin and Lomax [Bal] 
have been extensively tested, and shown to perform well in the incompressible, transonic 
and supersonic regimes [Vl,Yl], The Baldwin and Lomax model has the advantage that it 
does not require to calculate the boundary layer thickness, and has been widely used and 
tested. However, for detached or near detachment boundary layers, both the Cebeci-Smith 
and the Baldwin-Lomax model give erroneous results; this is due to the fact that they both 
rely on Prandtl’s mixing length theory, which is no longer relevant after detachment. In 
separated regions, we will use another algebraic model, introduced by Goldberg [Gol], 
which prescribes analytically the values of k and e in a separation bubble. 

All these models have been previously used and validated; the point of this paper 
is double: we first want to show that turbulent flows can be efficiently calculated on 
unstructured grids, using algebraic models; this is not straightforward because the direction 
normal to the wall is no longer a coordinate axis, so that nonlocal effects are not easy 
to compute; the details of the implementation are given in part II. Secondly, we want 
to investigate how well these algebraic models, which were defined after incompressible 
theories, and thus ignore density fluctuations, will perform for high speed flows. In part 
III, we will compute a boundary layer at Mach 7.4, and compare the results with the 
experiments of Hopkins et al. [Hoi]. No data on hypersonic flows over compression 
corners is known to the author; however, such a high speed computation part of the 
Hermes workshop [Mai], will be performed in part IV. 

II. ALGEBRAIC TURBULENCE MODELS AND UNSTRUCTURED GRIDS 

Linear algebraic turbulence models describe the Reynolds turbulent stress tensor a t as 

a t = n t D{u) (2.1) 

and the turbulent heat flux q t as 

1 , = ( 2 . 2 ) 

where Pr t is the turbulent Prandtl number, usually assumed constant, and where n t is the 
eddy viscosity, which is an explicit function of the flow variables W, although this function 
is usually not local. More precisely, 

o, Vo) = f{W ( x 0 , y) , y > 0) (2.3) 
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The value of /x t at a given point depends of the flow variables at all the points at the same 
streamwise location. 

One of the main advantages of algebraic models, which makes them so popular, is their 
minimum requirement of computer time and storage. Indeed, when one uses an t, j grid 
(Fig. 2) , the value of the eddy viscosity on the grid, /x,-,-, can be calculated cheaply since 

M«„i, = /(HW > 0) (2.4) 

and, in (2.4), many terms in Hu 0 jo are common to all points of same streamwise location, 
and so can be computed only once for each value of t<j. 

Whether these desirable features can be conserved for triangular grids, which can be 
of type a, but also of type b (Fig. 3), is the question that must be answered. 

A triangular grid can be, and usually is, at least locally, not orthogonal to the wall. In 
our approximation ([Rol]), the flow variables are linear on each triangle, and consequently, 
we will look for nodal values of the eddy viscosity. The normals to the wall cannot, as 
in the rectangular structured case, be approximated by a sequence of nodes. There are 
two ways to calculate the field of eddy viscosities: the first is, for each node, to draw the 
normal to the wall passing through this node, to track all the elements it gets through 
(Fig. 4), and then, using flow variables which are piecewise linear on this normal, compute 
the eddy viscosity, which is valid only for the considered node, since it is the only node 
with this streamwise location. This is very expensive, both in storage, since the list of the 
elements used for each node must be stored, unless it is recalculated at each time step, 
leading to an enormous cpu cost, and in cpu, since a complete viscosity calculation on a 
normal must be performed for each node. 

The other way, which we will use, is to consider a discrete set of normals, namely 
the normals to the wall, passing through the middle of the wall edges (Fig. 5), compute 
and store the elements they get through, together with the barycentric coordinates of 
the middle of their intersection with each element, and then, at each time step, use that 
information to compute the eddy viscosity on each normal. For each node, it is then 
necessary to interpolate the values of the eddy viscosity from the two neighboring normals 
(Fig. 5). 

More precisely, the path, in terms of elements and edges, followed by the normal, is 
computed through an algorithm first derived for characteristic methods ((Dal]), taking in 
account all the geometrical possibilities: the normal hitting a node, coinciding with an 
edge, etc. This is done once, at the beginning of the calculation; the result which consists, 
for each wall edge, of a list of triangles, and of the barycentric coordinates of the middle of 
the intersection of the normal with each triangle, is stored. The memory needed for that 
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is equal to the number of wall edges multiplied by the average length of a normal (which is 
fixed by an a priori bound to the boundary layer thickness) times one integer, the number 
of the element, and two reals, the coordinates. 

Stol = Wall Edges x Average Normal Length x (II + 2 R) (2.5) 

To be able to efficiently make the interpolation of the eddy viscosities from the normals 
to the whole mesh, one must also store some geometrical quantities. For each node in the 
boundary layer, one must store the number of the normal just on its left (say), its distance 
to this normal, plus one integer and one real to allow positioning on this normal, and one 
integral and one real to allow positioning on the normal on the right (Fig. 6) . 

Altogether, the memory needed to sore the interpolation data is 3 integers and 3 reals 
for each of the nodes in the boundary layer. 

Sto2 = B.L. Nodes x (3 1 + 3 R) (2.6) 

The global storage needed for the turbulence is 

Sto = Wall Edges x Average Normal Length (11+ 2J?) + Boundary Layer Nodes (31+ 3#) 

(2.7) 

For example, a flat plate calculation using 100 nodes on the wall and an average 25 nodes 
in the boundary layer will require 90 Kbytes of storage for the turbulence data, which is 
quite reasonable. The generation of this data takes about 40 seconds on a SUN workstation 
for this case. 

Once this data is generated, at each time step, the eddy viscosity is evaluated on the 
normals. This is now straightforward, since we have a normal represented by a set of 
points at which the flow properties, including the vorticity, are known; we are brought 
back to the structured orthogonal case. Then the interpolation is performed, which is also 
straightforward, since all the needed coefficients are stored. Altogether, the calculation 
of the eddy viscosity when one uses the Baldwin-Lomax model takes less than 1.5% of 
the global cpu time for an implicit scheme, and less than 4% for an explicit scheme. The 
evaluation of the eddy viscosity for the aforementioned case takes about 1.2 seconds on a 
SUN. 

The interpolation procedure described earlier is valid only when the solid body is 
convex; when its not, a given point of the flowfield is usually considered as influenced 
by two or more walls, and the influences are proportional to the inverse of the distance to 
the wall. A similar interpolation can be performed; the only added storage is: one real for 
each node of the boundary layer, representing the ratio of the influence of the eventual two 
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relevant walls, and in sto2, a boundary layer node must be counted twice if it is influenced 
by two convex components of the wall. This can be handled in a completely automatic 
and geometry independent way, which is consistent with the finite element spirit. 

Although the geometrical part will be much more complicated, the same thing can be 
done in three dimensions. Normals are drawn starting from the wall faces, and then a 
spatial interpolation is performed between these normals. 

in. HIGH SPEED ATTACHED BOUNDARY LAYERS 

Turbulent boundary layers have been studied quite extensively in the last two decades, 
and described in much details [Cel], Many authors have proposed turbulent closure models 
to predict these boundary layers; the simpler ones, which will be our concern, are those 
which suppose the turbulent stress tensor —pu'v 1 to be proportional to the deformation 
tensor D: 

—pu'v' = Ot = p t D(u) (3.1) 

where p t is the eddy viscosity, which in algebraic models is a function of the flowfield. 

Many eddy viscosity laws p t = Pt{W) have been proposed through the years, the 
successful ones separate the boundary layer in an inner and an outer part, which behave 
differently. In the inner part, the main feature of the analysis is Prandtl’s mixing length 
theory, which predicts the eddy viscosity to be 

pt = pl 2 \u\ (3.2) 

where |cj| is the magnitude of the vorticity and l the mixing length. According to Prandtl’s 
theory, in the fully turbulent regime, l is proportional to the distance to the wall: 

l = Ky (3.3) 

To account for the laminar sublayer, this expression has to be modified in the near wall 
region, as we will see later on. 

In the outer part of the boundary layer, it is generally admitted that the eddy viscosity 
is almost constant, although different values of this constant have been proposed. Cebeci 
and Smith [Cel] suppose that 

p t = ccu e 6i (3.4) 

where a is Clauser’s constant, u e the speed at the edge of the boundary layer, and 5, 
the incompressible displacement thickness. Again, this expression must be modified when 
approaching the edge of the boundary layer to account for the relaminarization of the flow. 
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From these ideas, Cebeci and Smith proposed an eddy viscosity model defined as: 


Ht = P& M 

£ = /cy(l-exp(-y/A)) ► 
for y < y c 


(3.4) 


and 


fit = pau^ii 


7 =(1 + 5.5 (y/5) 6 )- 1 \ 


(3.5) 


y>y e 


where 1 — exp (—y/A) 
the laminar sublayer, 
given by: 


is Van Driest’s damping factor, which allows the representation of 
A is the damping length, which according to Cebeci and Smith is 


A = A+ 



A + = 26. 


(3.6) 


7 is Klebanoff’s intermittency factor, which accounts for the alternatively turbulent and 
laminar regime which occurs around the boundary layer edge. It has been established 
experimentally by Klebanoff [KU]. The crossover value y e is obtained by requiring the 
continuity of fi t . 

This model has been tested extensively by its authors, and shown to perform well in 
the incompressible, subsonic and transonic case. Its main drawback is the necessity to first 
define, and then calculate, the boundary layer edge location, and the properties of the flow 
at that point, in order to calculate the outer eddy viscosity. This is very difficult to do 
numerically with a reasonable accuracy. 

To overcome this difficulty, Baldwin and Lomax [Bal] proposed an alternate outer 
formulation 


Pt = ocC ep pY MAX F MA x 7 

y>y e 


(3.7) 


where 

F{y) = y|w|(i ~ «P(- j)) (3.8) 


and Fmax is the maximum of F in a profile, and Y MAX the value of y at which it takes 
place. They take 


6 = 


Ymax 

Ckl 


(3.9) 


6 



where C ep and Ckl are constants determined by (3.9) and (3.10): 

C cp YmaxFmax = u e 6i 


(3.10) 


for theoretical profiles. York and Knight [Yl] have shown that, although the values C cp = 
1.2 and Cki = 0.65 could be used for low speed flows, for high speed flows, these values 
depended very strongly not only on the Mach number, but also on the skin friction, making 
the model quite unreliable. Another difficulty of the model is that for high speeds, the 
function F(y) does not exhibit a sharp peak as in the transonic case [Bal], making the 
determination of Ymax and Fmax difficult and unreliable. 

As remarked before, the difficulty involved in using the Cebeci-Smith model comes 
from the necessity to find u e , Si, and 6. We have: 

* = t\ 1 - -)dy 

Jo u t 

(3.H) 

u e Si = / (u e — u)dy. 

Jo 

It is clear that an inaccurate determination of 6 and u e = u e (6) will result in a very 
inaccurate value of the eddy viscosity. However, integrating (3.11) by parts, we have: 

uA= H v % dv 

which for an attached flow is equivalent to: 

f s 

u t 6i = / y\u\dy (3.12) 

Jo 

For high Reynolds number flows, the vorticity |w| always decreases sharply as the distance 
to the wall increases, even when the boundary layer experiences pressure gradients or 
interaction with a shock. The function y|u;| also decreases quite quickly, as we will see in 
the numerical results, so that to calculate u e 6, using (3.12), we only need a rought estimate 
of 6. We obtain this estimate through the variations of the Baldwin and Lomax function 
F: we stop the integration in (3.12) at the point y* where F(y*) is lower than Fmax * 
where (3 is an arbitrary constant. Since the function y|w| decreases quite quickly, (3 = 0.5 
gives an accurate enough result, and allows to separate between the boundary layer and 
an eventual shock. 

u e Si = j y|u;|dy. (3.13) 

To calculate the intermittency factor 7, we need a more accurate estimate of 8. We 
define a scaling length by 

_ Stv'Wv (3 14) 

v “~ irvMdy' 1 ' 
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By comparison with the theoretical profiles of Sun and Childs ([Sul]), it is found that: 


6 = YJC kl 


(3.15) 


where C K l is a constant which depends only slightly on the flow parameters, see Fig. 7, 
and which can be taken to be: 


Ckl = 0.45 


(3.16) 


for all practical purposes. 

In the original Baldwin and Lomax model, the damping length A was calculated using 
only wall quantities. 


A = A + 


P-vi 

yj PwT w 


There is no significant difference with (3.6) for low Mach numbers, because the density 
and molecular viscosity do not vary much, but for high speed flows, there is a difference, 
and Cebeci’s expression (3.6) was found to give better results. 

Finally, the model we are using for zero pressure gradient boundary layers is fi t given 
by (3.4) in the inner layer, with the damping length given by (3.6), and in the outer layer, 


Pt = pau e 6r y 


(3.17) 


with u e 6i given by (3.13), and 

1 = (1 + 5.5(2p^) e )-‘ (3.18) 

Vav 

where the length scale y av is given by (3.14). 

Two high speed flows over flat plates were computed, and the results were compared 
with the experimental data of Hopkins et al. [Hoi]. Unfortunately, high speed measure- 
ments are quite scarce, and only skin friction values were available. The first case is a flow 
at Moo — 7.4, Reoo/m = 8x 10 6 . The plate has a length of 2 and a height of 0.35; the grid 
has 2960 nodes and 5688 elements (Fig. 8), it was generated by splitting the quadrangles 
of a 74 x 40 grid through their diagonal. The shock is captured at the leading edge. The 
free stream temperature is Too = 97.3 K, while the temperature at the wall is T w = 311 A", 
which is about one third of the adiabatic temperature. The molecular viscosity is given 
by Sutherland’s law. 

Figure 9 shows the pressure coefficient, skin friction and heat flux versus the streamwise 
coordinate x. The capturing of the shock is done sharply, as one can see from the low 
deviation of the pressure coefficient (less than 1%) (the leading edge is at x = 0). The 
skin friction and heat flux have the awaited behavior, starting at very high values and 
decreasing sharply as the Reynolds number Re x = Re w • x increases. The temperature of 
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the wall is lower than the adiabatic temperature, so that the heat is transferred from the 
flow to the wall. On Figure 10, the skin friction is plotted versus the momentum Reynolds 
number, and compared with the data of Hopkins et al. The agreement is very good, except 
for the first point, at Re g = 800, which is in transitional flow, while no attempt to account 
for the transition has been made here. The streamwise variations of 0, u e x <5, and Y av are 
also presented. It can be seen that they are very smooth, although a few discontinuities 
appear near the end of the plate, where the higher and coarser part of the grid plays a 
role. The crosswind profiles at x = 1.45 of the speed, the density, the turbulent and the 
total shear stress, and of u> * y are presented on Figure 11. Although no experimental data 
was available for comparison, the results are quite reasonable, at least qualitatively. The 
profile of density exhibits a minimum shortly of the wall; this is consistent with the fact 
that the wall is colder than the corresponding adiabatic wall. This peak of density induces 
a peak in total shear stress, because the eddy viscosity is proportional to the density. This 
behavior is maybe not physical, and a drawback of the mixing length theory. We have 
checked that it is not dependent on whether local or wall flow properties are used in the 
calculation of the damping length. As announced before, the function u*y decreases quite 
quickly when approaching the boundary layer edge, and it can be integrated easily with a 
reasonable accuracy. This computation took about 10 hours of cpu on a Gould, of which 
only about 3% was due to the evaluation of the eddy viscosity. The main over cost when 
compared to laminar calculations comes from the slower convergence to steady state: the 
resolution algorithm is a linearly implicit pseudo-time marching to steady state, in which 
the variations of the eddy viscosity with respect to time are not taken in account: at time 
level n, the viscosity is frozen, and the algorithm is advanced of one time step, giving 
the flow properties at time n + 1, at which the eddy viscosity is re-evaluated. This still 
allowed the use of a Courant number of 50 without encountering stability problems, but 
the number of iterations necessary to converge to steady state was about 50% more than 
in the laminar case, giving a 50% overcost. 

Another of the test cases of Hopkins et al., at a higher wall temperature, was computed. 
The geometry is the same as previously; the Mach number is M 0 0 = 7.4, the Reynolds num- 
ber is Reoo/m = 3 x 10 7 , the free stream temperature is Too = 38K, the wall temperature is 
T w = 305 K. The grid used has 1170 nodes and 2204 elements; it was obtained by splitting 
in triangles a 30 X 39 grid. In Figures 12, 13, and 14, the same results as for the previous 
case are presented. Because of the much coarser streamwise discretization, the results are 
less smooth, and the agreement with experimental data is not as good as for the previous 
case. Nevertheless, the error on the skin friction (about 8% at the maximum) is less than 
the estimated experimental error, so that the agreement can be considered reasonable. 
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IV. ALGEBRAIC MODELS FOR HIGH SPEED SEPARATED FLOWS 


It is well-known that, when the boundary layer becomes separated, both the Cebeci- 
Smith and the Baldwin-Lomax model give erroneous results, because they rely on Prandtl’s 
mixing length theory, which is no longer relevant. This has been observed by many authors 
(see e.g., [VI], [Yl]), and is maybe the major drawback of these models. Some improvement 
can be obtained by using the Cebeci-Smith model, and modifying the Van Driest damping 
factor ([VI]), but the results are still not very good in and downstream of the separated 
zone. 

In 1986, Goldberg ([Gol]) proposed a new algebraic k—e model to account for separated 
regions. His model is based on the following assumptions and observations on the separated 
region: the stress scale is given by the maximum shear stress in the separated layer, not 
by the wall stress; the shear layer has qualitatively the same turbulent structure when 
it is detached as when it is attached, and the length scale is the height of the separated 
region. Considering this, and continuity arguments, he proposed to take the kinetic energy 
of turbulence to be 

1 — ex 
k = kf, 

where the subscript b refers to the backflow edge (defined by the point where the tangential 
speed is zero), and where <f> is a parameter, found empirically to be <f> = 0.5. 

The turbulence energy dissipation is taken to be 


p(-«K£) 2 ) 


Vb' 

JL\ 


(4.1) 


er- — 

yt 

because the length scale is y b . 

For high Reynolds number flows, k b can be taken to be, by analogy with the attached 
case, 

. u* 

ki = v& (4 - 3) 

where C* = 0.09, and where 

u t = V^-(uV) max 

if — (tt'v')mo* 18 the maximum turbulent shear stress, which occurs in the detached layer, 
and must be provided by the model used outside the separation bubble. 

Altogether, the kinematic eddy viscosity is taken as 



*4 = 

= Ci u $yb\[*ff{y) 


(4.4) 
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where 


ttv) = A?- + B ( 4 . 5 ) 

yt 

The function / accounts for the laminar part in the vicinity of the wall, the constant 
A and B are found to be optimal at 



c; = 0 . 7 . 


The function G is the Gaussian: 

«■«> 

Goldberg tested his model by computing a supersonic flow over a compression corner, 
and obtained impressively good results, including the correct prediction of separation and 
reattachment [Golj. 

To take in account separation, we will use the following blend of Cebeci’s and Gold- 
berg’s models. In the attached regions, we use the Cebeci model as described before, only 
replacing in the van Driest damping factor the wall shear stress by the maximum shear 
stress in the profile, to avoid an unphysical reduction of the eddy viscosity near separation 
and reattachment. 

For separated profiles, consistently with Goldberg’s hypothesis, we suppose that the 
shear layer is not qualitatively disturbed by the separation, so that we can take the eddy 
viscosity in it to be given by the Cebeci model, provided we use the distance to the backflow 
edge instead of the distance to the wall. In other words, n t is given by 


for yt<y <y e and 


for y >y c where 


Ht = pP\w\ 

t = x[y — y t ) (l - exp (-^)) 




( 4 . 7 ) 


( 4 . 8 ) 


( 4 . 9 ) 
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(4.10) 


_ jT y 2 \ u \ d y 
Vav /» yl w l d y ’ 

In the separation bubble, i.e., for y < yj, we use Goldberg’s model, as described by (4.4). 

The implementation of this model is no major problem in our framework, since we have 
well defined normals to the wall, on which we can easily compute yt, y*, and the integrals 
(4.9) and (4.10) giving u e 6,- and y av . The cost, in terms of cpu or storage, is not significantly 
different for this model as compared with the preceding one, and remains very small (« 
2% of the cpu time is used to compute the eddy viscosity) . 

Separation usually occurs at points where the solid wall is not convex; at these points, 
the influences of the two convex components on one given fluid point are taken into account, 


and averaged according to 


^2^1 + dj/X2 


(4.11) 


" d. + li, ' ' 

where the subscripts 1 and 2 refer to the two convex components of the wall, and where d 
is the normal distance to the wall. 

The program automatically recognizes convexity defaults, and makes the needed aver- 
ages. The parts of the eddy viscosity depending on the different convex components are 
computed separately, to allow vector processing; all this is completely geometry indepen- 
dent, assuming the solid wall is locally convex. 


As a preliminary test of this model, a high speed, moderate Reynolds number flow 
overa a 15° compression corner was performed. The grid is shown in Fig. 15, it has 
2782 nodes and 5345 elements, it was obtained by refining a cartesian grid, using an 
algorithm defined by Pouletty and Palmerio ([Pol], [Pal]). The shock is captured at the 
leading edge, the corner is at 1.39 meters from the leading edge. The Reynolds number is 
Reoo/m = 4.95 x 10 5 , the Mach number is Moo = 5; the free stream temperature is Too 
= 83. 6K, the temperature of the wall is Too = 288K. Fig. 16 is a plot of the pressure. It 
can be seen that both the bow shock and the corner shock are captured neatly. On Figure 
17, the pressure, skin friction and heat flux coefficients are plotted versus the streamwise 
coordinate (the corner is at x = 0). The boundary layer experiences a small separation 
in the vicinity of the corner, after which the skin friction quickly recovers high values, as 
awaited. Fig. 18 is a plot of the ratio n t / n of the eddy viscosity to the molecular viscosity; 
it varies between 0 and 15 between the leading edge and the corner; the flow in this region 
is first laminar and then transitional. After the corner nt/ n increases quickly to values 
around 100, indicating a fully turbulent flow. Fig. 19, 20, and 21 are plots of the velocity 
vectors at the corner, after separation and near the outflow respectively. Fig. 22 shows 
profiles of speed, density, and total and turbulent shear stress at the corner (x = 0). These 
profiles have the awaited shape, although the discretization is maybe a bit coarse. 
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V. CONCLUSION 


We have shown that algebraic turbulence models can be used in conjunction with 
unstructured grids, at no major overcost, both in terms of cpu time and storage; the 
program remains completely geometry independent, which is consistent with the spirit of 
finite elements and unstructured grids. 

A practical way to use the Cebeci-Smith model has been proposed, for both attached 
and separated flows; in the latter case, Goldberg’s modification has been used in the 
separated regions. The model has been shown to give accurate skin friction for high 
speed no pressure gradient boundary layers. A preliminary result for a separated flow is 
presented. 

The major remaining open problem to be solved before the model can be trusted to give 
even coarse results on real configurations is transition. Whether an ad hoc representation 
by switching the model of if the predicted value of the eddy viscosity is lower than a critical 
value, as suggested by Baldwin and Lomax [Bal], would give a relevant result, has not 
been investigated. More experimental data on transitional and/or separated high speed 
flows is certainly necessary before answers can be given. 


13 


References 


[Bal] B. S. Baldwin and H. Lomax, “Thin layer approximation and Algebraic model for 
separated turbulent flows,” AIAA 78-257, Huntsville, January 16-18, 1978. 

[Ca] G. V. Candler and R. W. MacCormack, “Hypersonic flow past 3D configurations,” 
AIAA 87-0480, Reno, January 12-16, 1987. 

[Cel] T. Cebeci and A. Smith, “Analysis of turbulent boundary layers,” Academic Press, 
1974. 

[Chi] S. Chakravarthy, “High resolution upwind formulations for the Navier-Stokes equa- 
tions,” Von Karman Institute Lecture Series, 1988-05, March 7-11, 1988. 

[Dal] F. el Dabaghi, Thesis, Universite Paris XIII, 1984. 

[Gol] U. C. Goldberg, “Separated flow treatment with a new turbulence model,” AIAA 
J., Vol. 24, No. 10, October 1986, pp. 1711-1713. 

[Hoi] E. J. Hopkins, E. R. Keener, T. E. Poleh, and H. A. Dwyer, “Hypersonic turbulent 
skin-friction and boundary layer profiles on nonadiabatic flat plates,” AIAA J., Vol. 
10, No. 1, January 1972. 

[Jl] A. Jameson and H. Rieger, “Solution of steady three dimensional compressible 
Euler and Navier-Stokes equations by an implicit LU scheme,” AIAA 88-0619, 
Reno, 1988. 

[Kll] P. S. Klebanoff, “Characteristics of turbulence in a boundary layer with zero pres- 
sure gradient,” NACA TN 3178 (1954). 

[Mai] M. Mallet, J. Periaux, P. Perrier, and B. Stoufflet, “Flow modelization and com- 
putational methodologies for the aerothermal design of hypersonic vehicles: Ap- 
plication to the European Hermes,” AIAA 88-2628, June 1988, San Antonio. 

[Osl] S. Osher and S. Chakravarthy, “Upwind difference schemes for the hyperbolic 
systems of conservation laws,” Math. Comp., April 1982. 

[Pal] B. Palmerio, “Self adaptive FEM algorithms for the Euler equations,” INRIA 
Report 338, 1985. 

[Pol] C. Pouletty, “Generation et optimisation de Maillages en elements finis, applica- 
tion a la resolution des equations de la mecanique des fluides,” Th&se de Docteur 
Ingenieur, Ecole Centrale, December 1985. 


14 


[Rol] P. Rostand and B. Stoufflet, “TVD schemes to compute compressible viscous flows 
on unstructured meshes,” Proceedings of the 2nd International Conference on 
Hyperbolic Problems, Aachen (FRG), 1988 (Vieweg). 

[Shi] J. Shang and W. L. Hankey, “Numerical solution of supersonic turbulent flow over 
a compression ramp,” AIAA J., Vol. 13, October 1975, pp. 1368-1374. 

[Spl] C. G. Speziale, “On nonlinear K -t and K-e models of turbulence,” J. Fluid Mech., 
Vol. 178, 1987, pp. 459-475. 

[Stl] B. Stoufflet and L. Fezoui, “A class of implicit upwind schemes for Euler simula- 
tions with unstructured meshes,” to appear in J. Comp. Phys. 

[St2] B. Stoufflet, J. Periaux, L. Fezoui, and A. Dervieux, “Numerical simulations of 
3D hypersonic Euler flows around space vehicles using adapted finite elements,” 
AIAA 87-0560, Reno, January 1987. 

[Sul] Sun, C. C. and M. E. Childs, “A modified wall wake velocity profile for turbulent 
compressible boundary layers,” J. Aircraft, Vol. 10, June 1973, pp. 381-383. 

[VI] M. Visbal and D. Knight, “The Baldwin-Lomax turbulence model for two- 
dimensional shock-wave/boundary layer interactions,” AIAA J., Vol. 22, No. 7, 
July 1984. 

[Yl] B. York and D. Knight, “Calculation of a class of two dimensional turbulent bound- 
ary layer flows using the Baldwin-Lomax model,” AIAA 85-0126, Reno, January 
1985. 


15 



16 



17 




19 




20 



Fig. 6. 
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Fig. 7. Variations of Ckl with the flow parameters. 
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